Phylogenetic analysis and genetic diversity of the xylariaceous ascomycete Biscogniauxia mediterranea from cork oak forests in different bioclimates

Cork oak is a tree species with ecological importance that contributes to economic and social development in the Mediterranean region. Cork oak decline is a major concern for forest sustainability and has negative impacts on cork oak growth and production. This event has been increasingly reported in the last decades and seems to be related with climate changes. Biscogniauxia mediterranea is an endophytic fungus of healthy cork oak trees that turns into a pathogen in trees weaken by environmental stress. Understanding the drivers of B. mediterranea populations diversity and differentiation is expected to allow a better control of cork oak decline and preserve forest sustainability. Endophyte isolates from different cork oak forests were identified as B. mediterranea and their genetic diversity was evaluated using phylogenetic and microsatellite-primed PCR analyses. Genetic diversity and variability of this fungus was correlated with environmental/phytosanitary conditions present in forests/trees from which isolates were collected. High genetic diversity and variability was found in B. mediterranea populations obtained from different forests, suggesting some degree of isolation by distance. Bioclimate was the most significant effect that explained the genetic variability of B. mediterranea, rather than precipitation or temperature intensities alone or disease symptoms. These findings bring new implications for the changing climate to cork oak forests sustainability, cork production and quality.

Phylogenetic trees of each single locus were not able to individually resolve B. mediterranea isolates for any factors of interest (forest, bioclimate, disease severity index or disease symptoms; results not shown). The concatenated dataset comprised sequences from 52 B. mediterranea isolates and from the three outgroup species  19,20 , but no evident patterns regarding factors of interest were exposed using this approach (Fig. 1). Isolates obtained from cork oaks growing in different locations and distinct bioclimates did not clustered together; and visible tree disease symptoms did not contribute for the distribution of isolates. For example, Clade A comprised isolates (7; e.g. Bm25, Bm57, Bm37) from different forests and bioclimates (from the most humid to the driest), also presenting isolates obtained from trees with different disease severity levels. The same was observed for Clade C (with 6 isolates; e.g. Bm36, Bm66, Bm08), which even included an isolate (Bm79) obtained from an olive tree. Moreover, some information can be inferred from the tree. Group III (with moderate ML support) displayed a higher number of isolates from trees with mild symptoms (21 out of the 31 found in Group III, contrasting with 7 isolates in the remaining 19). Clade B (in Group I) only included isolates (Bm50, Bm54 and Bm69) from declining trees, while Clade D (Group III) only included isolates (Bm55, Bm11 and Bm03) collected from mild diseased trees, though both clades presented isolates from different forests and bioclimates. Overall, isolates from different locations were distributed along the phylogenetic tree, with few of them being placed together. The same distribution pattern was observed regarding isolates from different bioclimates, even though the sequences in Group VII refer to isolates (Bm49, Bm47, Bm46 and Bm64) obtained from two forests (GV and GR) with the same bioclimate (sub-humid). Also, a clear distribution pattern was not observed concerning the presence of exudates on the B. mediterranea host trees. However, subclade A1 only encompassed isolates (Bm25, Bm23, Bm57 and Bm37) obtained from cork oak trees not producing trunk exudates. The opposite was observed for Clade B (Bm50, Bm54 and Bm69) and Group VI (Bm55, Bm11, Bm03 and Bm15), where all isolates were recovered from trees producing exudates, even if from dissimilar bioclimatic locations (sub-humid and humid). In line with the low resolution of B. mediterranea isolates by any of the studied factors (forest, bioclimate, disease severity index or disease symptoms), the majority of B. mediterranea isolates collected from the same tree were not clustered together. For instance, although pairwise identity within those isolates was high, isolates from a single tree, such as Bm60, Bm61, Bm63 and Bm67 (pairwise identity of 95.4%) and Bm46, Bm56 and Bm58 (99.6%), were found to cluster better with isolates obtained from other forests and/or were placed apart in the phylogenetic tree. Similarly, isolates collected on olive trees (Bm79 and Bm80) were distantly placed from each other, in separate lineages.
Microsatellite-primed PCR fingerprinting of 68 B. mediterranea isolates generated different banding patterns. Primer (GTG) 5 generated 15 bands (from 0.25 to 1.5 kb), (CAG) 5 generated 17 bands (ranging from 0.3 to 2 kb, in which one was monomorphic), (ACAC) 5 generated 17 bands (from 0.25 to 1 kb) and M13 generated 26 bands (from 0.15 to 1.5 kb). The monomorphic band was removed and a binary dataset with the remaining (74) band positions was used for molecular analysis. Using this approach, the genetic diversity of B. mediterranea varied within different factors of interest/populations (Table 1). Considering total populations, the number of alleles (Na) varied from 1.40 to 1.97 and the number of effective alleles (Ne) from 1.45 to 1.57, which were found considering 'forest' and 'exudates' factors, respectively. These results were corroborated by the genetic diversity found within populations (Hs; Table 2). 'Exudates' factor (total population) also displayed the highest Shannon's information index (I = 0.50), Nei's gene diversity (h = 0.33) and percentage of polymorphic loci (PPL = 98.7), while 'forest' factor shown the lowest genetic diversity (I = 0.38; h = 0.26; PPL = 67.9). In addition to 'exudates' , other disease-related factors also revealed high genetic diversity levels for total populations, in particular when considering 'cankers' and 'disease severity levels' factors (Tables 1 and 2 www.nature.com/scientificreports/ (Table 1). Considering the 'bioclimate' factor, isolates obtained from sub-humid forests revealed the highest genetic diversity and those from hyper-humid the lowest. While Grândola (GR, a sub-humid forest) was the forest with the highest B. mediterranea genetic diversity, PG-ER (hyper-humid) and AL (humid) were the forests with the lowest diversity. Isolates collected from healthy trees (low disease severity level) revealed the lowest genetic diversity and declining trees were associated with more diverse B. mediterranea isolates. This result is in line with the highest genetic diversity found among B. mediterranea isolates collected from trees producing trunk exudates. However, the opposite was found for trees with trunk cankers, in which higher genetic diversity  www.nature.com/scientificreports/ was associated with isolates from trees without this symptom. Regarding defoliation, the isolates obtained from trees with very accentuated damage were also less genetically diverse, while those from trees with moderate and light damages presented higher genetic diversity.
In agreement with phylogenetic tree analysis from multilocus sequence analysis, principal coordinates analysis (PCoA) did not revealed clusters of cork oak isolates, according to the factors of interest (Fig. S1, both axes only capturing around 15% of data variation). Also, isolates from olive tree did not cluster all together and were closer to isolates from healthy and mild diseased trees. However, certain populations revealed a higher differentiation in allele frequencies than others, as evaluated by Gst coefficient (Table 2). There was a higher genetic differentiation among 'forest' (Gst = 0.239), followed by 'defoliation' (Gst = 0.151), and 'bioclimate' (Gst = 0.097) populations. The lowest genetic differentiation was found in 'cankers' (Gst = 0.046), 'disease severity levels' (Gst = 0.045) and 'exudates' (Gst = 0.032) populations. As expected, gene flow (Nm) values were opposite to Gst, being higher among populations with lower genetic differentiation and vice-versa (Table 2). These results agree with the genetic pairwise distances found in such populations (Table S1). Pairwise genetic distances in 'forest' achieved higher values (ranged from 0.244 to 0.067), than 'defoliation' (0.200 to 0.019), or bioclimate (0.085 to 0.053). Other disease-related variables ('cankers' , 'disease severity levels' and 'exudates') never achieved more than 0.048 genetic distances. Among populations, isolates obtained from sub-humid bioclimate were the most genetically close from the ones obtained in the other forests (Table S1). Indeed, GR (a sub-humid forest) displayed the least genetic distance from all other forests. In contrast, isolates from humid bioclimate (in particular from AL forest) displayed the highest genetic distance from all other forests. Interestingly, isolates obtained from very accentuated defoliated trees revealed a high genetic distance from other isolates. AMOVA results revealed that variation within populations was always higher than among populations, which reinforces the high variability of B. mediterranea composition. In any case, variation among populations was higher for 'forest location' (10%), followed by 'bioclimate' (6%), 'cankers' (4%) and 'exudates' (3%), all at p < 0.001 (Table 3), suggesting higher genetic differentiation between regions and bioclimates. Accordingly, Mantel test shown significant correlation between B. mediterranea genetic diversity and geographic distance of cork oak forests (R = 0.105, p = 0.001). For further understanding the relative contribution of forest location, environmental, and disease-related factors in driving the genetic diversity of B. mediterranea populations, a redundancy analysis was performed (Fig. 2). The variables that explained the variation of B. mediterranea genetic diversity were 'forest location' and 'bioclimate' factors, followed by 'exudates' and 'cankers' (all at p < 0.001), being all the others (temperature, precipitation, disease severity levels, and defoliation) not significant. The combination of all significant variables explained 9.602% (p = 0.001) of B. mediterranea genetic variance. Most of this variation is due to the 'bioclimate' and 'forest The production of exudates showed higher correlation with 'bioclimate' (0.188%) than with 'forest location' (0.081%). The opposite was found for the presence of cankers, in which no correlation was found with 'bioclimate' , but a 0.943% correlation was found for 'forest location' variable. The index of association for the estimation of linkage disequilibrium revealed that LI and GV forests, as well as forests from the humid, sub-humid and semi-arid bioclimates was significantly different from 0, which indicates linkage disequilibrium (Table S2).

Discussion
In this study, the genetic diversity of endophytic B. mediterranea isolates was evaluated by multilocus phylogeny and microsatellite fingerprinting (MSP-PCR). The evaluated isolates comprised those obtained from cork oak trees, thriving in different forests, bioclimates, and displaying distinct charcoal disease severity levels and symptoms. This experimental approach is expected to provide new information about which factors/variables contribute the most to this species diversity and variability. Other isolates obtained from olive trees were also analyzed to evaluate host-specificity. As expected, all studied B. mediterranea isolates were phylogenetically apart from other Biscogniauxia species included in this study as outgroups, revealing the genetic divergence of B. mediterranea species. In general, individual (i.e., single-locus) and concatenated (i.e., multi-locus) phylogenetic analyses did not correlate B. mediterranea isolates with any of the studied environmental/disease variables. Also, after determining the microsatellite polymorphic patterns, B. mediterranea isolates did not clustered according to the factors of interest. These results suggest a high genetic diversity and variability of B. mediterranea populations, which is concordant with other reports 9,15,17,18 , even when isolates come from the same stroma 16 . Furthermore, the inability to resolve B. mediterranea isolates obtained from cork oak and those obtained from olive tree reinforces current knowledge about the high adaptability of this fungus to different hosts, as suggested by 17 .
The genetic variability of B. mediterranea isolates was mostly explained by 'forest location' (6.9%) and 'bioclimate' (4.6%). Indeed, when these factors were combined, they explained 8.8% of B. mediterranea genetic variation. Considering 'forest location' , B. mediterranea populations demonstrated the lowest genetic diversity within population, while also revealing the highest genetic differentiation and lowest gene flow among populations. These results suggest some isolation by distance of B. mediterranea communities, which agrees with the significant correlation (Mantel test) found between the B. mediterranea genetic diversity and forest geographic location of cork oak forests from which they were obtained. This contrasts with previous studies, where B. mediterranea intraspecific polymorphism and genetic diversity were not associated with geographic position of host trees 17 . However, the contribution of geographic isolation in B. mediterranea genetic differentiation was suggested for cork oak populations in Tunisia 18 . The high genetic variability of B. mediterranea has been related with the heterothallic mating system of this species 15 and sexual reproduction with the production of a high number of variable ascospores 21 . Our results suggest that global B. mediterranea population displays linkage disequilibrium and a clonal genetic structure. Despite that, populations from the majority of forests (PG-ER, PG-RC, AL, GR, HC-CT and HC-MA) demonstrated to have random mating with frequent sexual reproduction. B. mediterranea ascospores are primarily dispersed by wind after the occurrence of precipitation 21,22 , although insects could also play a role for their spreading in short-and long-distances, depending on their bioecology 21,22 . Our results suggest low migration rates between geographically distant B. mediterranea populations, indicating short-distance dispersal. The clonal structure of the populations from the majority of bioclimates and the finding that PG-ER www.nature.com/scientificreports/ and AL forests presented the lowest genetic diversity among all other sampled forests reinforces this suggestion.
In contrast with other sampled forests, both are characterized by a high density of mixed forest trees and low anthropogenic interference, which may restrain spore dispersal to and from distant locations. Indeed, canopy architecture and the use of mixed tree species have been reported as a management strategy to reduce spore dispersal of pathogenic fungi [23][24][25] . In addition, high genetic differentiation among forest populations can be a result of random events but there is enough gene flow to refute the effects of genetic drift 26 . While a significant variation was found among forest populations, high variability of B. mediterranea was also described within populations. This result concurs with the high intraspecific genetic variability found within populations described by Henriques et al. 27 , indicating an adaptation of B. mediterranea to the environment and ensuring species longterm survival.
Besides 'forest location' factor, a significant variation on B. mediterranea populations occurred in response to 'bioclimate' , although precipitation and temperature alone were not significantly correlated with B. mediterranea variability, as revealed by variation partitioning redundancy analysis. In contrast, other studies indicated a positive correlation between B. mediterranea genetic diversity with temperature and rainfall 18 . The ability of B. mediterranea to develop in a wide range of temperatures 27 , associated with the significant variability of this fungus in different bioclimates may represent a problem for charcoal cork oak disease management. This will be further challenged by the effect of combined bioclimate and charcoal symptoms (exudates in cork oak trunk) in increasing B. mediterranea variability. Indeed, in Tunisia, the correlation between bioclimate and charcoal disease development has been suggested 7 . In any case, the variability promoted by disease symptoms (exudates and cankers presence) was better explained when taking into consideration the forest location, suggesting that other factors characteristic of each forest (i.e., silvicultural practices not included in this study) are contributing to B. mediterranea variability. Moreover, B. mediterranea isolates collected from declining trees or trees with charcoal symptoms, like trunk exudates, presented a higher genetic diversity than those collected from healthier trees.

Conclusions
Several reports suggested that Quercus suber populations and their associated microbiota are vulnerable to different bioclimates [28][29][30] and will be affected by the predicted climate changes 31 . Therefore, as cork oak forests currently displaying a moist bioclimate become more arid, they will be increasingly affected by environmental stressors. Taking all together, the results obtained in this work support the previous suggestions that B. mediterranea isolates have facilitated adaptation. Our findings reinforce the previous knowledge of B. mediterranea opportunistic behavior and reveal the importance of bioclimate as a source of B. mediterranea variability, exacerbating the implications of a changing climate on cork oak forests sustainability which will affect cork production and quality.

Methods
Biscogniauxia mediterranea isolates. Biscogniauxia mediterranea isolates were obtained from cork oak twigs as endophytes 32 . Twigs were sampled, during 2017, from eight Portuguese cork oak forests with different locations and bioclimate classifications ( Fig. S2; Table S3). Bioclimate for each forest was determined using Emberger index and Emberger climatogram 33,34 and ranged from hyper-humid to semi-arid bioclimates [hyper-humid (PG-RC and PG-ER), humid (AL and LI), sub-humid (GV and GR) and semi-arid (HC-MA and HC-CT)]. If existent, trees with different disease severity levels were collected from each location. Evaluated disease symptoms included defoliation (5 levels: 0-10%-no damage; 11-25%-light damage; 26-50%-moderate damage; 51-90%-severe damage; > 90%-extreme damage), as well as canopy and trunk damages (3 levels: 0no damage; 1-moderate damage; 2-severe damage) (Table S4). Disease damages included dried, wilting and decolorated leaves, presence of cankers, decolorated trunk and presence of exudates (Fig. 3). Disease severity levels were grouped into three categories (healthy, mild, and declining), considering the combination of different symptoms and corresponding levels, as described elsewhere 32 . Endophytic isolates were obtained by sterilizing the surface of cork oak twigs, plating twigs onto Potato Dextrose Agar (PDA) medium and obtaining pure cultures through re-plating outgrowing fungi into fresh PDA medium 30 . DNA of pure cultures was extracted using Quick-DNA Fungal/Bacterial Miniprep Kit (Zymo Research, Irvine, CA, USA) and B. mediterranea isolates were identified using universal primer pairs ITS1-F (5′-CTT GGT CAT TTA GAG GAA GTAA-3′) and ITS4 (5′-TCC TCC GCT TAT TGA TAT GC-3′) 35 . A total of 74 B. mediterranea isolates obtained from cork oak were distinguished and identified by sequencing (Table 4). Three other isolates were obtained from Olea europaea twigs (Bm78 36 ) and olives (Bm79 and Bm80 37 ), collected from cvs. Cobrançosa (Bm78 and Bm80) and Madural (Bm79). These O. europaea-derived isolates were included in this study to evaluate host-specificity in B. mediterranea. All methods complied with relevant institutional, national, and international guidelines and legislation.
DNA sequences were processed using AB1 trace files in Geneious version 2010.4.8.5 (Biomatters, New Zealand), unless stated otherwise. For each isolate and molecular marker, forward and reverse sequences were trimmed (0.05 error probability limit), assembled and consensus sequences were created. Consensus sequences obtained in this study were deposited in GenBank (for accession numbers, see www.nature.com/scientificreports/ alignments of each region were performed by using MUSCLE version 3.5 algorithm with a maximum of 10 iterations. Distance measure used for 1st iteration was kmer4_6 and subsequent iterations were run with pctid_ kimura. All iterations were performed using UPGMB clustering method and CLUSTALW sequence weighting scheme. If needed, alignments were manually edited and Gblocks (web-based, version 0.91b, January 2002, http:// molev ol. cmima. csic. es/ castr esana/ Gbloc ks_ server. html, last accessed date: 08/02/2021) was used to eliminate poorly aligned positions and divergent regions, allowing smaller final blocks 38 . Geneious version 2010.4.8.5 was then used to concatenate alignments. Some isolates were not sequenced with enough quality in some targeted loci and were excluded from individual alignments before concatenation. DNA sequences from closely related taxa-Biscogniauxia atropunctata, B. nummularia and Xylaria hypoxylon 39,40 -were used as outgroups and were retrieved from GenBank (Table 4). Phylogenetic trees for each individual locus were generated using sequences from 74 B. mediterranea isolates for ITS, 70 for TEF, 68 for GS, 66 for ACT , 69 for CHS and 71 for TUB2 ( Table 4). The final concatenated alignment used to build the multilocus phylogenetic tree included sequences from 52 B. mediterranea isolates and from the three outgroup species (isolates marked with * in Table 4). PartitionFinder2 version 2.1.1 was run on CIPRES Science Gateway (web-based, version 3.3, https:// www. phylo. org/ porta l2/, last accessed date: 17/02/2021) 41 to find best-fit partition schemes of each loci 42 . Bayesian inference (BI) trees were computed using Markov chain Monte Carlo (MCMC) algorithm in MrBayes version 3.2.7 43 . Models (lset) and prior probability distributions (prset) were set according to PartitionFinder2 results. Two independent runs were performed with one million generations and four chains in each run. The two runs were converged with a burnin of 25% and tree with posterior probabilities was generated. Maximum Likelihood (ML) trees were generated using W-IQ-TREE (http:// iqtree. cibiv. univie. ac. at/, last accessed date: 17/02/2021) 44 , a web server for IQ-TREE 45 . Best-fit model was computed using ModelFinder version 1.4.2 46 , with an edge-linked partition model 47 . Branch support analysis was performed using 1000 ultrafast bootstrap replicates 48  Microsatellite-primed PCR fingerprinting. The molecular diversity of 68 B. mediterranea isolates was also evaluated by using microsatellite-primed PCR (MSP-PCR) with four sets of primers, (GTG) 5 , (CAG) 5 , (ACAC) 5 and M13 (phage M13 core sequence; 5′-GAG GGT GGNGGNTCT) 17,18 . PCR reactions were prepared   Table 4. GenBank accession numbers for each sequenced locus [internal transcribed spacer (ITS), translation elongation factor 1-α (TEF), partial glutamine synthetase (GS), actin (ACT ), chitin synthase 1 (CHS) and β-tubulin 2 (TUB2)] of endophytic B. mediterranea isolates obtained from cork oak twigs and other fungi used in this study. Information is given concerning the sampled forests and disease severity level of sampled host trees. Isolates without GenBank accession number were only used for microsatellite-primed PCR fingerprinting. a Refers to sequences used for multilocus analysis. b Refers to B. mediterranea specimens obtained from olive trees.   53 . Given the reduced sample size, B. mediterranea isolates collected from olive tree (Bm78, Bm79 and Bm80) were included only in phylogenetic and PCoA analyses. Redundancy analysis (RDA) was used to explore responses of B. mediterranea composition to environmental (bioclimate, mean maximum and minimum temperatures and mean total precipitation for the 10 years previous to sampling collection) and disease variables (disease severity level, exudates, cankers, and defoliation), by making use of the R version 4.0.2 54 . Analyses were performed using the package vegan version 2.5-7 55 , except when stated otherwise. Spatial trend was included in RDA using a trend surface analysis. Latitude-longitude data was transformed into flat Cartesian coordinates using geoXY() of SoDA package (version 1.0-6.1) 56 . To compute polynomials of degree 3, poly() of STATS package (version 4.0.2) was used. Spatial and environmental variables were analyzed separately on a first approach. RDA was performed using rda(), while anova.cca() was used to perform Monte Carlo permutation test (1000 permutations) and test significance of global model. Forward selection of explanatory variables was performed using ordistep() with 999 permutations, using spatial and environmental variables to find most parsimonious RDA. Then, another RDA was performed, as described before, for selected variables. The contribution to variation of those variables was performed using RsquareAdj() and Monte Carlo permutation test (1000 permutations) to determine significance. For explaining the variation on B. mediterranea composition, variation partitioning was performed for the best model, using the most explanatory variables (combination of spatial and environmental variables), and making use of varpart(). Partial-RDA was performed to evaluate the influence of conditional variables obtained from variation partition. Statistical significance was performed as described for RDA. Spatial variable was referred as 'forest location' variable in sections related with this analysis. Linkage disequilibrium of B. mediterranea populations was evaluated by the index of association (I A ) and the unbiased index of association ( r d ), in which clonal populations have results significantly different from 0 and sexual populations do not have statistical significant I A and r d 61 . Analysis was performed using popsub() to specify population and ia() with 999 permutations for calculation using poppr package (version 2.9.0) 62 .
Received: 1 August 2021; Accepted: 31 December 2021 Table 5. Locus regions amplified for B. mediterranea phylogenetic analyses and correspondent PCR conditions. The corresponding amplicons size are shown in brackets.